<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of NM_example</title>
  <meta name="keywords" content="NM_example">
  <meta name="description" content="">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html v1.5 &copy; 2003-2005 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../index.html">Home</a> &gt;  <a href="index.html">CO4OI</a> &gt; NM_example.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for CO4OI&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->

<h1>NM_example
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong></strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>This is a script file. </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre class="comment"></pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
This function is called by:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
<!-- crossreference -->



<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre>0001 clear;
0002 
0003 
0004 addpath(genpath(<span class="string">'/tptool'</span>)) <span class="comment">%change this path to the Tensor Product Toolbox path in your computer</span>
0005 addpath(<span class="string">'operators'</span>)
0006 addpath(<span class="string">'Data'</span>)
0007 addpath(genpath(<span class="string">'misc'</span>));
0008 addpath(genpath(<span class="string">'Solvers'</span>));
0009 
0010 
0011 <span class="comment">%% Generate image example</span>
0012 
0013 <span class="comment">% input parameters</span>
0014 n= 8; <span class="comment">% n rows, image dimension n x m</span>
0015 m= 8; <span class="comment">% m columns</span>
0016 N=n*m;
0017 plots=1; <span class="comment">%auxiliary flag for plots</span>
0018 noisy=1; <span class="comment">%auxiliary flag indicating presence of noise. 1: noisy case; 0: noisless case</span>
0019 input_snr=30; <span class="comment">%input snr</span>
0020 
0021 <span class="comment">% parameters image</span>
0022 <span class="comment">% field type:'image' loads image from file (natural image resized) // 'spikes' creates a synthetic image with k spykes</span>
0023 paramIM.dim=[n m];
0024 paramIM.type=<span class="string">'spikes'</span>; 
0025 paramIM.k=6;
0026 paramIM.F0=1;
0027 paramIM.fileimage=<span class="string">'eta-carinae_16x16.mat'</span>;
0028 
0029 <span class="comment">% create an image in x</span>
0030 x=gen_astro_object(paramIM); 
0031 <span class="keyword">if</span>(plots), figure, imagesc(x),colorbar, axis image, <span class="keyword">end</span>
0032 xu=x(:); <span class="comment">%unfolded signal</span>
0033 Tx=fast_tensor_product(xu); <span class="comment">% outer product of the image vector -&gt; Tx is a (n*m, n*m, n*m) tensor</span>
0034 
0035   
0036 <span class="comment">%% Generate mask for the open frequences</span>
0037 u=1;
0038 mask=nonredundant_mask(N);
0039 
0040 sigmau=1/4;sigmav=1/4;sigmaw=1/4;
0041 T=mask_3D_3(n,m,u,sigmau, sigmav, sigmaw);
0042 Tr=build_redundant_table( T );
0043 
0044 maskB=zeros(N,N,N);
0045 <span class="keyword">for</span> i=1:size(Tr,1)
0046   I=Tr(i,:);
0047   maskB(I(1),I(2),I(3))=1;
0048 <span class="keyword">end</span>
0049 
0050 
0051 
0052 <span class="comment">%% Add noise only to the non-redundant part of the measurements vector</span>
0053 
0054 y0=analop_fft_6d(Tx,n,m); <span class="comment">%compute 3-dimension 2D-Fourier Transform</span>
0055 yB=y0(maskB&amp;mask); <span class="comment">%measured values, only once</span>
0056 NB=numel(yB);
0057 
0058 
0059 
0060 idmask=find(maskB&amp;mask==1);
0061 
0062 <span class="keyword">if</span> noisy
0063     <span class="comment">% Add Gaussian i.i.d. noise</span>
0064     
0065     eB=sqrt(1/NB*sum(abs(yB(2:end)).^2));
0066     param.sigma_noise=.5*10^(-input_snr/20)*eB;
0067     
0068     nBr=(randn(size(yB))); 
0069     nBi=(randn(size(yB)));
0070     
0071     nBr=sqrt(NB)*param.sigma_noise*nBr./norm(nBr(:),2);
0072     nBi=sqrt(NB)*param.sigma_noise*nBi./norm(nBi(:),2);
0073     
0074     nB=nBr + 1i*nBi;
0075     yB=yB+nB;
0076     
0077     <span class="comment">% noise bound based on a chi-squared model</span>
0078     param.epsilon=sqrt(6)*sqrt(2*NB + 4*sqrt(NB))*param.sigma_noise ;
0079     <span class="comment">% tolerance on the noise bound</span>
0080     param.epsilon_up= sqrt(6)*(sqrt(param.sigma_noise^2*(2*NB) + 5*param.sigma_noise^2*sqrt(2*(2*NB))));<span class="comment">%expectation + 2.1 std</span>
0081     param.epsilon_low= sqrt(6)*(sqrt(param.sigma_noise^2*(2*NB)  -5*param.sigma_noise^2*sqrt(2*(2*NB))));
0082 <span class="keyword">else</span>
0083     param.epsilon=1e-7; <span class="comment">%default 1e-7 or less, depending on the scale of the signal (ideally should be 0)</span>
0084 <span class="keyword">end</span>
0085 
0086 <span class="comment">% build measurement vector</span>
0087 maskB=logical(maskB);
0088 Ty=build_full_tensor(yB,maskB&amp;mask); <span class="comment">% symmetrization of the tensor (replicate measurements)</span>
0089 y=Ty(maskB);
0090 y=y(:);
0091 
0092   
0093 <span class="comment">%% Parameters of the algorithm</span>
0094 
0095 <span class="comment">% general parameters</span>
0096 param.max_iter = 200;
0097 param.rel_obj=5e-3;
0098 param.verbose = 2;
0099 param.N=N;
0100 
0101 <span class="comment">% specfic parameters L2-ball projection</span>
0102 param.A = @(x)analop3dfft(x,n,m, maskB); 
0103 param.At = @(x)synop3dfft(x,n,m, maskB);
0104 param.nuB2=1;
0105 param.tightB2=1;
0106 param.verboseB2=2;
0107  
0108 <span class="comment">% specific parameters NN-prox</span>
0109 param.C=@(x)C(x,N);
0110 param.Ct=@(x)Ct(x,N);
0111 param.nu_NM=N;
0112 param.tight_NM=1;
0113 param.max_iter_NM=200;
0114 param.rel_obj_NM=5e-3;
0115 param.verbose_NM=2;
0116 
0117 <span class="comment">%estimate initial point</span>
0118 Aty=param.At(y);
0119 param.xinit = Aty;
0120 <span class="comment">%compute threshold parameter for the NN-prox</span>
0121 aux=sv(Aty,N); <span class="comment">%docu!</span>
0122 param.lambda_NM=0.1*max(abs(aux))/param.nu_NM; 
0123 
0124 
0125 
0126 <span class="comment">%% Solve NM problem</span>
0127 t_start=tic;
0128 
0129 sol=solve_nm(y,param);
0130 
0131 <span class="comment">% projection (real)positive orthant</span>
0132 sol=real(sol); sol(sol&lt;0)=0;
0133 
0134 <span class="comment">% rank 1 approximation</span>
0135 [ x_NM] = s_hopm( sol, N, 1e-3 );
0136 x_NM=x_NM./sum(x_NM(:)); <span class="comment">%normalisation</span>
0137 x_NM=reshape(x_NM,n,m);
0138 
0139 t_end=toc(t_start);
0140 
0141 <span class="comment">%% Log</span>
0142 
0143 maxsnr=calculate_snr(x,x_NM); 
0144 
0145 <span class="comment">% plot the result using the same scale as for the original plot</span>
0146 mx=max(x(:));
0147 <span class="keyword">if</span> (plots), figure, imagesc(x_NM, [0,mx]),colorbar, axis image,  <span class="keyword">end</span> 
0148 
0149 fprintf(<span class="string">'\n M/N= %e, SNR0 = %e, Ellapsed time= %e\n\n'</span>, u, maxsnr, t_end);
0150 
0151 
0152 
0153 
0154 
0155</pre></div>
<hr><address>Generated on Thu 18-Jul-2013 19:25:15 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>